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The Kohn variational method of calculatina scatterina phase shifts 
has been examined oreviously, and calculations of the phase shifts of 
S-wave electrons on atomic hydrogen done, by Schwartz. This paper 
oresents the results of using a trial function somewhat simpler than 
Schwartz's and compares the two sets of results. In addition this 
paper presents and compares the results of two different normalizations 
of the variational principle and gives a somewhat more detailed 
development of the variational principle used than in most previous 



papers. 
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I. INTRODUCTION 



In 1961 Schwartz used Kohn's variational principle to calculate 
the scattering ohase shift of S-wave (4 = 0) electrons incident on 
atomic hydrogen. The trial wave function he used had the form 



with as many as 50 terms in the sum [Ref. 1]. 

This paper oresents the results of calculations for the same 
problem using a different trial wave function. Assuming the results 
of Schwartz's calculations are correct, a comparison of the new results 
to his is made to determine the relative merits of the new trial wave 
function. The new calculations were also done using two different 
normalizations of the variational principle, and a comparison between 
these results is also presented. As a lead-in to the actual calcula- 
tions this paper first presents a development of the theory behind 
them. 







where 
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II. THEORY 



A. VARIATIONAL THEORY 
1 . General Form 

First in the presentation of the theory is a general look at 
the variational method of determining the exact wave function of a 
system. The general variational principle starts with the integral 



T- dT which equals 0 when 4^ is the 

exact wave function. After some manipulation it is found (See Appendix 
A-1), in spherical coordinates, that 

SI= 2^ JSY(£-H)YdT I 

- llm f r’-nir YJ-n- 

or SI = lim [ f - J 



when is the exact solution of h 4^= Ef. 

2. Specific Forms 

Now, by specifying the radial portion of 4^ at r = 0 and as 
r-^oo the general result above yields a particular variational prin- 
ciple. For example, the necessary conditions imposed on T for a bound 
state of angular momentum i are that at r = 0 ?'^(r) = 0 where = 

3od that as r-^oo y-^O. These give O 0 which 
is the Rayleigh-Ritz variational principle for bound states. The 
Kohn variational method for scattering phase shifts is derived from 
this general formula by putting in a set of conditions appropriate to 
the scattering problem. Therefore an examination of the appropriate 
conditions is now presented. 
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B. SCATTERING THEORY 



1 . Phase Shifts and Partial Have Analysis 

First note that the phase shift being calculated is the phase 
difference between the radial portions of the asymptotic forms of the 
wave functions with the scattering potential and without the scattering 
potential. Dividing the incoming wave into partial waves by angular 
momentum (partial wave analysis) there exists a phase shift for each 
partial wave. This set of phase shifts completely determines the 
scattering cross section. In the case of low energy scattering all 
the phase shifts, except the S or = 0 wave phase shift, are very 
small, and therefore the S-wave phase shift is often the only one 
needed for accurate calculations. 

2. Asymptotic Conditions 

So the conditions needed are those on the radial portion of 
the solution for angular momentum i of (E - H)4'= 0 in spherical 

/A. 

coordinates where H is the Hamiltonian with the scattering potential 
included. At r 0 the radial portion must be finite 

which implies that must equal 0. As r-^°o , and the potential 
approaches must become like the wave function of a free particle, 
and the most general form ^ can assume in spherical coordinates is 
shown by Schiff [Ref. 2] to be 

•h h 

where is the phase shift of the ^ partial wave, and is 

an arbitrary amplitude function. 
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C, FINAL VARIATIONAL PRINCIPLE 






1 . Final Form and Normalization 

Putting these conditions into the general form and allowing 
to be the variable parameter for yields the specific form 



= ~ "f ^ (See Appendix A-1). 

Now choose fC"^^) = secant(-^jj) and write "J' ^ 

so that ~ - ^ j^gi^ 

S(l ■*■ P('3^^)/k ) - ^ (l which is Kohn's variational 



principle. 



Ordinarily Kohn's principle is written [tan(7^^)/0 = J +• 
which means that when 'f is exact = tah('^7a)/k So by 



t« 



nl2^ 



taking the variation of I + with respect to several variable 

parameters in a trial function 4^, one of which and finding the 

values of these parameters which make I i" stationary, 

+ tariC'^^)/]^) = 0, one gets first order values of and the 
other oarameters. These first order values are then substituted back 
into 1+ to give the stationary result L^^'n('^f)/k] . 

2. Usefulness of Variational Theory 

What good is all this? It means that to find the scattering 
phase shifts all one has to do is select a reasonable trial wave 
function with several variable parameters which satisfies the above 
conditions, put it in the integrals, find parameter values that satisfy 
S(l+t«7*n('%)/k) = 0 and put these values back into [ta^(V^)/l(I = 

X+tctnf'%)/k If one can find a relatively simple and easy-to-integrate 



( 

wave function that gives good results for a difficult problem, this 
method may save considerable work. This is a big if with complications. 
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III. CALCULATIONS 



A. GENERAL OUTLINE 

In the case of variational calculations of scattering phase shifts 

A 

a unique problem appears. In scattering calculations E - H has a 
continuum of eiqenvalues which passes through 0. In the calculations 
E - H is represented by a matrix M-. in the subspace of the complete 

* J 

Hilbert space of H which is spanned by the N terms in the trial function 
This matrix has N eigenvalues and when any of these eigenvalues becomes 
very small tan(7j|^)/h may become very large and inaccurate. Thus at a 
given energy E convergence to better values with more terms is not 
guaranteed. However, by adding a non-linear variable K to the trial 
wave function and varying it for each value of E the behavior of 

can be mapped out. Then the best value of tan(->jf^)/K is 
found by picking the value from the flattest part of the graph of 
tan( 9 ^^)/(^ versus /C. As the trial wave function becomes more and more 
accurate with more and more variable parameters the flat regions of 
the curve become flatter and longer; converging, hopefully, to better 
and better values of the phase shift. When the improvement becomes 
small compared to the size of tan('^^)/K the final value is determined 
from the flat region of the graph and the probable further improvement 
due to more terms. The error is estimated from reading the graph and 
how well the further improvement can be determined. A good presenta- 
tion of this general method of calculation is given in Schwartz [Ref. 3] 
The method used in the present case differs by using the values from 
the curves directly without estimating probable improvement and 
estimating the error only from reading the graph. 
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B. THE HAMILTONIAN 



Now apply the above to scattering S-wave electrons off atomic 
hydrogen. In this problem there is an incoming electron, the hydrogen 
proton and the electron bound to the proton, all of which interact with 
one another. To take in all these interactions the following 
Hamiltonian was used 






-2^ - 



where radial distance and energy are in Bohr radii and Rydberg units 

, ■2. 

respectively (See Appendix B). Therefore R is the energy of the 
incoming electron, and -1 is the energy of the bound electron in the 
ground state. To avoid the complications of excitation of the bound 
electron k was limited to less than .866, which corresponds to about 
10.2 ev, so only low energy scattering was considered. 



C. THE TRIAL WAVE FUNCTION 

However one still had to take into account the indistinguishability 
of the electrons, and this was done with the trial wave function. The 
following wave function was used 












/\ 

where is the parity operator that interchanges r.| and and where 



^0 



ikr^ ) = 






. • *1 , Smkir; Coskr, 

and j, (kr,) = are spher- 



ical Bessel functions. The -f sign gives a symmetric wave function which 
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implies spin angular momentum S = 0, and the - sign gives an anti- 
symmetric wave function which implies S = 1. Thus the 1 ~ P-|2 
care of the indistinguishability of the electrons and their fermion 
properties . 

is the asymptotic part of the wave function and approaches 



and 



#1 jecar.tC'K) e"'"' 

^fTT ( ir^ 



as r,— >«»«» 



as r 2 ~^<>o 



contains the necessary asymptotic plane wave condition for one 
electron, coupled with a ground state condition for the other electron 
in the e"*^ term. The sum portion contains the close-in terms with the 
variable parameters Note that is also a variable parameter. 

D. INTEGRALS AND FINAL FORMULAS 

Now with two electrons two radial coordinates, r.| and r 2 » and 
their difference, r.| 2 » are used. Does this change the variational 
principle? The answer is no, with the proper normalization of the 
trial wave function (See Appendix A-2). So the variational principle 
remai ns = I + . Now defining 



Ri= r° + :»rI 

Mu= JTJTx 
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and 



one solves ^Ct 



VN)- iMuCi = 



0 



5 ( I ^/l^) ~0 j;;o gg-f- ^f^g gp^j ^ ^ ordeT. Theo 

inserting these values back into the variational equation gives 
_ . - _ ( W, t B, t 1/K ) 

L - Wo + Bo ~ + B,.) 

where Wg = c9r9, W, = 2^c]r 9, and = ^c]r] as shown by 

Schwartz [Ref. 3] . These are the final integrals and equations used 
in the calculation of the phase shift using f(‘^) s secant(->j>) . The 
phase shift was also calculated using f{n^) ~ cosecant(-)^) which gives 

fcot(’)^)/K] , J - cot(yj/k and requires 



(I ~ e 

The final formula is 

- ( W B, “ l/fe)_ 

L VkJ ^ Wk+Bj.' ^(Wo^-Bo) 

where - cot('>^) and the W's and B's remain the same as before except 
B - + -h B 2 , C.^;\c9fcl, and R^.= TiR® + r]. 

E. THE ACTUAL CALCULATIONS 

The actual calculations are done as follows. First the integrals 

Bq, B-| , B 2 , r9, r| , are done by hand for the general case, and then 

the final numbers are obtained by computer evaluation of the general 

integrals for particular cases. Then Tm..C. = -R. is solved independ- 

j ij J 1 

ently of the first order solution for /> since C. - c9 c| and M. . 

independent of ^ implies yM^-^-c9 = -r 9, and .cl = -rL This 

T ' J J f K 

solution is most accurately done by matrix inversion. The 
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factors in the close-in terms of 4^ are there to increase the accuracy 
of the matrix inversion and have no effect on the value of the phase 
shift. Finally the W's and the stationary value of-^are calculated 
using the previously given formulas. 

The above calculations were done for values of fc from 0.6 to 2.5 
in steps of 0.1 for each value of k from 0.1 to 0.8 in steps of 0.1 
with up to 40 close-in terms. This was done for both the symmetric 
and antisymmetric wave functions and both secant and cosecant normal- 
izations. Then the values of tan('>j>)/K or -cot(->^)/K were plotted 
versus for each case by the computer (See Figs. 1, 2, 3, 4), and 
values of tan(‘»^)/t< and -cot{n^^)A\ selected from the flattest portion 
of each curve. The actual values of *>^ were then calculated from these 
numbers for comparison with Schwartz's values. Additional values were 
obtained for 0.4 at 49, 56, and 72 terms to see how much better 
the values might get and where significant inaccuracies in the matrix 
inversion would occur. 
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Fig. 1 - Graph of tan(Y)/K versus K (kappa) fork=0.5 and S = 1. 



18 



TAN»f/k 




Fig, 2 - Graph of tan(*^)/K versus if (kappa) for K= 0.5 and S 
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-COT»j/k 




Fig. 3 - Graph of -cot(-^)/K versus ic (kappa) fork= 0.5 and S = 1 . 
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Fig. 4 - Graph of -cot(‘>^)/|^ versus K (kappa) for k= 0.5 and S=0. 
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IV. DISCUSSION OF RESULTS 



A. COMPARISON WITH SCHWARTZ'S RESULTS 

Table I lists the values of tan('>j?)/K and Yobtained from the 
present calculations and from Schwartz's calculations for S = 1. Table 
II does the same for S - 0 and Tables III and IV do the same for 
-cot(Y)/k for S “ 1 and S * 0 respectively. The uncertainties in the 

present values listed in the tables represent the errors one would 

estimate from the graphs in the absence of Schwartz's values for compar- 
ison. Comparison of the results is done quite easily from the tables, 
and the following results are obtained. In Table I the percentage 
differences between the values of present calculations for^ and 
Schwartz's calculations are between 2% and 8 % with one case of 12%. 

In Table II the differences in ^ are between 10% and 20% with one 
case of 40% and one of 207%. In Table III the differences in ^ are 
between 2% and 8% with one case of 15%, and in Table IV the differ- 
ences are between 9% and 20% with one case of 30% and one of 213%. In 

the cases with 207% and 213% differences in^ the differences in the 

absolute values are only 2% and 9% respectively but the values have 
opposite signs. These two cases are the only such cases of difference 
in sign, both occur at K- 0.3 for S = 0, and neither can be explained 
adequately at the present time. Ignoring these two very large cases, 
the average differences from Schwartz's values in ^ both secant 
and cosecant normalizations, are 5.5% for S = 1 and 15.9% for S ~ 0. 
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Table I. List of present secant normalization values versus Schwartz's 
values for S = 1. The numbers in parentheses give the uncertainty in 
the last digit. 





Present 


Schwartz 


K 


TanC-^)/!^ 


7 ^( radians) 


Tan('?p)/K 


')j>( radians) 


0.1 


-2.460(5) 


-0.212(5) 


-2.056(4) 


-0.2028(4) 


0.2 


-2.555(5) 


-0.472(1) 


-2.260(3) 


-0.4245(5) 


0.3 


-2.745(5) 


-0.6880(9) 


-2.492(4) 


-0.6419(8) 


0.4 


-3.125(5) 


-0.8960(8) 


-2.833(2) 


-0.84776(35) 


0.5 


-3.723(8) 


-1.0778(9) 


-3.384(3) 


-1.0370(4) 


0.6 


-5.85(6) 


-1.293(3) 


-4.40(1) 


-1.2087(7) 


0.7 


-8.35(2) 


-1.4013(4) 


-6.74(2) 


-1.3619(6) 


0.8 


-28.75(5) 


-1.52734(8) 


-17.2(6) 


-1.498(3) 



Table II. List of present secant normalization values versus Schwartz's 
values for S = 0. The numbers in parentheses give the uncertainty in 
the last digit. 





Present 


Schwartz 


k 


Tan (i^)/t^ 


>^(radians ) 


Tan(7^)/K 


^>p( radians) 


0.1 


-10.00(5) 


-0.785(2) 


-6.68(2) 


-0.589(1) 


0.2 


-16.5(3) 


-1.276(5) 


-9.23(2) 


-1.0743(9) 


0.3 


31.5(5) 


1.465(2) 


-26.4(1) 


-1.4452(5) 


0.4 


8.15(5) 


1.273(2) 


15.87(4) 


1.4145(4) 


0.5 


3.59(1) 


1.062(1) 


5.17(2) 


1 .2017(13) 


0.6 


2.135(1) 


0.908(1) 


2.845(8) 


1.041(1) 


0.7 


1.455(1) 


0.7946(3) 


1.917(5) 


0.9303(12) 


0.8 


1.10(2) 


0.721(9) 


1.530(5) 


0.8857(16) 
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Table III. List of present cosecant normalization values versus 
Schwartz's values for S = 1. The numbers in parentheses give the 
uncertainty in the last digit. 





Present 


Schwartz 


k 


-Cot(^)/M 


o^( radians) 


Tan(77)/K 


•>y( radians) 


0.1 


41.920(5) 


-0.23417(3) 


-2.056(4) 


-0.2028(4) 


0.2 


10.10(3) 


-0.4597(12) 


-2.260(3) 


-0.4245(5) 


0.3 


4.1201(1) 


-0.68023(1) 


-2.492(4) 


-0.6419(8) 


0.4 


2.0505(2) 


-0.883859(5) 


-2.833(2) 


-0.84776(35) 


0.5 


1.09847(3) 


-1.068541(11) 


-3.384(3) 


-1 .0370(4) 


0.6 


0.572425(5) 


-1 .2399641(27) 


-4.40(1) 


-1 .2087(7) 


0.7 


0.25708(6) 


-1 .39275(4) 


-6.74(2) 


-1.3619(6) 


0.8 


0.056993(7) 


-1.525233(6) 


-17.2(6) 


-1 .498(3) 



Table IV. List of present cosecant normalization values versus 
Schwartz's values for S = 0. The numbers in parentheses give the 
uncertainty in the last digit. 





Present 


Schwartz 


k 


-Cot ( 77 ) /k 


T^(radians) 


Tan('>^)/h 


•Y( radians) 


0.1 


11.245(3) 


-0.72686(13) 


- 6 . 68 ( 2 ) 


-0.589(1) 


0.2 


1 .6855(2) 


-1.24566(4) 


-9.23(2) 


-1.0743(9) 


0.3 


-0.1188(1) 


1.53517(3) 


-26.4(1) 


-1.4452(5) 


0.4 


-0.77845(2) 


1.268932(7) 


15.87(4) 


1.4145(4) 


0.5 


-1.1076(2) 


1 .06504(7) 


5.17(2) 


1.2017(13) 


0.6 


-1.29527(3) 


0.910137(11) 


2.845(8) 


1.041(1) 


0.7 


-1.39933(3) 


0.79574(1) 


1.917(5) 


0.9303(12) 


0.8 


-1.409835(5) 


0.725378(2) 


1.530(5) 


0.8857(16) 



24 



B. COMPARISON OF DIFFERENT NORMALIZATIONS 

The average difference between the values of calculated with 
secant and cosecant normalizations is 2%. One other observation is 
worthy of note, and that is that the values of -cot(y)/k, found by the 
cosecant normalization have more apparent significant figures than the 
values of tan('>^)/i^ . This is due to the fact that the plots of -cot('^)/K 
versus /c do not show fluctuations in -cot(^)/^ anywhere as large as 
in tan(Y)/k • To see this compare Fig. 1 with Fig. 3. Since the 
normalization of the variational principle is arbitrary there should 
be no difference in accuracy from one to the other, and the comparison 
with Schwartz's results shows this to be true. Without Schwartz's 
results for comparison one would seem to have more confidence in the 
results of the cosecant normalization due to the larger number of 
apparent significant figures. However, without a standard of compari- 
son, there is no reason why one normalization could not give more 
significant figures in a given case and still not give any better an 
answer for 

C. ACCURACY AND CONVERGENCE 

The results examined above are for 40 terms in the sum and they show 
fairly good results for S = 1 and generally poorer results for S - 0. 

An examination of runs at 0.4 with 49, 56, and 72 terms shows the 
following. Beyond 40 convergence of the flat region toward Schwartz's 
values is very slow and not always guaranteed, since for the S = 0 
■secant normalization at 49 the value actually became worse. In addition, 
at 56 and 72 terms, significant errors in the matrix inversion were 
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noticed dna the accuracy beyond about 50 terms became questionable. 

The limit imposed by computer accuracy was therefore between 50 and 
60 terms, and within this limit continued convergence did not give 
appreciably better results. 

D. RELATIVE ACCURACY OF S = 1 AND S = 0 

The triplet (S = 1) calculation is better behaved than the singlet 
(S = 0) calculation, probably for the following reason: the anti- 

symmetric (triplet) 4^ is automatically 0 for r^ = r^ so that the 
electrons cannot coalesce. This effect is qualitatively the same as 
the effect of electrostatic repulsion, and thereby operates to remove 
some of the burden from the variational machinery. The singlet wave 
function contains no such built-in convenience. Thus for S = 1 the 
calculations have a greater absolute accuracy than for S = 0 and so 
should have a greater relative accuracy for S = 1 than for S ■= 0. 
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V. CONCLUSIONS 



There are two differences between the trial functions used in this 
paper and those used by Schwartz. In the asymptotic part Schwartz 
used instead of and in the close- 

in terms Schwartz included an r ^2 which was left out in the present 
calculations to make the integrals easier. From the comparison of the 
present results with Schwartz's results one concludes that in general 
the convergence of the curves for the present trial function is not 
as good in an equivalent number of terms. The idea that the effect 

p 

on convergence of eliminating the r ^2 factor could be offset by using 
more terms of a simpler form seems to be cut short by the inaccuracies 
in the matrix inversion and slower rates of convergence with larger 
numbers of terms. Therefore one concludes that a more complex trial 
function must be used to get good results, and that one can be misled 
by the apparently good accuracy (the small uncertainties in the tables) 
of the simpler trial functions without comparison. Finally one can 
conclude from the above comparison that the simpler form works 
reasonably well for S = 1 but not so well for S = 0, and that there 
is really no difference in accuracy between the two normalizations 
despite the larger number of apparent significant figures for the 
cosecant normalization. 
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APPENDIX A 



THE VARIATIONAL PRINCIPLE 



This appendix shows the steps in the derivation of the variational 
principle used, for the case of one radial coordinate r, and for the 
case of two radial coordinates r-j and 
1 . One Radial Coordinate 

Start with I = f H'i £-”H) f JT where E is a constant 

' s/ 

energy term and H = ' ^ U(*^) 

where U(r) is a central potential. Then write 

S1-- + Xt(E-H)ir <it] 

Now us1nq f 5 7' ( S' VS<(') “ V- ( jfV')') + 

one gets ^ 2 ^ Irlrf E-W^JT + JT 

- J'V'fJtV't'j Jr. 

u = 2^1 jH'fE-Hirjr + Tj 



= 2 ^ fiTCE-fiiv'jr 

since the integrals are over all space and dS- dJl- . 

si= liM 1..V, 

Ir-^ec 

for 4^ exact. With properly normalized the integral over the angular 
part can be made unity. Then 

\r^oO 

where Rjj is the radial portion of Ifl Then inserting the condition 
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exact. Note that in 



on at co 

_ . rr X ^in U^r - ^ 

r-^oo r 

one finds, using S R< “ R( 

that S2- “ f ^je)/fe for 

the case of one radial coordinate the correct normalization of for 

2. Two Radial Coordinates 

In this case H = ( V| +Vx) f M ('^ i ) 



Then as before 

r|-)oC> 

+ 1 i 1 3^. s 't'diiiJir - 1 1 wX.Xs,':" f jAidi: 

forM^ exact. Again with'f' properly normalized the integrals over 
and can be made unity so 

n = I i>n r f ^ ^ dr Rs] •! + 1 ^ " *1 ^ Ri 

l^-^CTO ° ^->0O 

Then inserting the conditions 

P — ’ ■ff-n'l e'"’- p -> rfi, I 

R«r,^«= Tt^«) kr, ^ J n‘>jj 

one gets ^ I = — ^ f where the|~can be eliminated 

by putting an additional factor of.^in the normalization of . So 
finally SI- “ f (yt)lk as before. Note that the 

correct normal ization of 4^ for 0 is 



now 
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APPENDIX B 



THE HAMILTONIAN 



IJrite HT-- t i * 4) 



then [-^ +sM 7f --t •* ^ 



)Y 



where 



- € 



is the ground state energy of hydrogen. Therefore 



[s-SJ 4 - = [ H- v,^ V." ^ 44-' ( 



2^ - J_ 



ri ' 4 






MS 



Then let -j^i - >^1 j j^x~ where <=*n- so when 

or one Bohr Radius. So 



f. ‘ft ' 1 



r= 



£k'^- 4?- + ( i " A ' 

t it'L JiLl" 

and dividing through by and setting K - 

[ - I f- i- ( fi ^ ^ 4^-0 

k 1^1 g'* ^ _ I ^ n 

Now E - TTJT* ' 2. “ K i\o where is the Rydberg energy 

, 2 . 

so I? measures energy in Rydberg units. Finally 

and for excitation to the first excited state E “ required 



which implies 



^ and 



k<S - . 



866, 
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CCMPUTER PROGRAM OF CALCULATIONS 



IMPLICIT PFAL*8(A-H,0-Z) 

INTEGER S( 100,2) ,T(100,2) 

REALMS K(10),KAPPA{25),MTX(100, 100) ,MTXI ( 100 , 100 ), LAMB 
1DAI2) 

DIMENSION R1 ( 100) ,R2( ICO) ,C1 ( 100 ) ,C2( 100) ,N( 5) ,P( 2), PH 
lASEI 2) 

NIT = 5 
NOT = 6 

READ(NIT,AOO) NSFT 
400 FCRMAT(I4) 

PEAD(NIT,3P9) (M I ) ,I = 1,NSET) 

399 F0PMAT(I4) 

PEADINI T.4C1 ) NP,NK,NKAPPA 
:*01 FCPMAT(3I4) 

READ! NIT. 402 ) (K(I ),I=1,NK) 

402 F0RMAT(Fl0.4) 

PE AD (MIT, 402) (KAPPA! I ) ,I=1,NKAPPA) 

READ(NIT,^C2) ( P ( I ) , I =1 , NP ) 

N1=N(NSET) 

2 READ(NIT,102) ( S( I ,1),T(I,1),S(I,2),T(I,2),I=1,N1) 

102 FrRMAT(4l4) 

DO 40 L1=1,NP 
DO ^0 L2=1,NK 

C CALCULATION OF INTEGRALS 

CALL PI (K(L2) ,B1,P(L1) ) 

CALL BII (K(L2) ,B2,P(L1) ) 

CALL BIII(K(L2) ,B3,P(LD) 

DC 40 M=1,NSET 
N1=N(M) 

WPITE(NOT,300) N1 , P ( L 1 ) , K ( L2 ) 

300 FORMAT! • 1 •, ‘CALCULATIGNS FOR*,' N = *,I4,4X,*P = »,F6. 
12,4X, *K = • ,F10.4/) 

WRITE! NOT, 100) 

100 FORMAT! //T26, • I • , 6X , • S ( I ) • ,4X , • T ! I ) • / ) 

WPITF!N0T,200) ! I , S ! I , L 1 ) , T ! I , L 1 ) , I =l , N1 ) 

20C F0RMAT!T25,I4,4X,I4,4X,I4) 

WPITE!M0T,103) 

103 FORMAT! //T24, 'KAPPA* ,14X, ' TAN/K * , 1 5X , * PHASE' , 14X, '-COT 
1/K* ,15X, 'PHASE * /) 

DO ^0 L3=1,NKAPPA 
DO 3 1=1, Ni 

3 CALL PI!$!I ,L1) ,T!I,L1),K!L2),KAPPA!L3),R1!I ),R2! I),P! 
ILD) 

DC 6 I=1,N1 
DO 6 J=1,I 

CALL MIJ(S!I,L1),T!I,L1),S!J,L1 ) , T ! J , L 1 ) , K ! L 2 ) , KAPPA! L 
13),MTX!I,J),P!LD) 

IF!J.EO.I) GO TO 6 
MIX! J,I )=MTX! I , J) 

6 CONTINUE 

C MATRIX INVERSION AND SOLUTION OF M*C = -R 
CALL GAUS31! Nl ,C. C,MTX,MTXI ,KER,100) 

DC 13 I=1,N1 
SUM1=0. 0 
SUM2=0.C 
DO 12 J=1,N1 
TEPM1=MTXI(I ,J)’»R1! J) 

TERM2=MTXI ! I ,J)*R2! J) 

SUM1=SUM1+TERM1 
12 SUM2=SUM2+T£RM2 
Cl !I )=-l«SUMl 
•13 C2! I )=-l*SUM2 
C CALCULATION OF W'S 
SUM1=0.0 
DO 14 1=1, Nl 
TERM=C1( I )*R1! I ) 

14 SUM1=SUM1+TERM 
W1=SUM1 
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SUM1=0.0 
DO 15 1=1, Ml 
TFRM=C 2 ( n*Rl( I ) 

15 SUM1=SUM1+TERM 
W2=2*SUM1 
SUM1=C.0 

or 16 1=1, Ml 
TEFM=C2( I )*R2( I ) 

16 SUM1=SUM1+TERM 
W3=SUM1 

C SECANT NORMALIZATION 

LAMBDA! 1 ) = W1 + B 1 - ( W2 + B2 1 1 /K ( L2 ) ) **2/ ( A * ( W3 ♦ B3 ) » 
PHASE! 1)=DATAN(K( L2)*LAMBDA( 1 » ) 

C COSECANT MOPMALI ZATIDN 

LAMBDA! 2)=W3+B3-! W2 + B2-1/K! L2M**2 / !A*! Wl + fil ) ) 
PhASF!2)=DATAM!-l / !K!L2 )* LAMBDA !2) )) 

WP ITE! MOT,201 ) KAPPA !L3) ,! LAMBDA! I ) » PHASE! I) t I = l» 2) 
201 F0RMAT!T20,F10.A,A!4X,D16.8n 
AO CONTI HUE 
END 



FL'NCTIOM FACTRL!M) 
IMPLICIT PEAL*8! A-H,0-Z) 
IF!M.LE.l) GO TO 2 
N = M-1 
PPOD = M 
DC 1 1=1, N 

1 PEOD = PRnD*!M-I) 

FACTRL = PPOD 

PE TURN 

2 FACTRL = 1.0 
RETURN 

END 



FUNCTION SIGF!M,A,J,K,B,C) 
IMPLICIT PEALA8! A-H,C-Z) 

PPOD = J 

IF! J.EO, !K + in GO TO 2 
I = 1 

1 PPOD = PROD*!J-I) 

IF! ! J-I ).EC. !K+1) ) GO TO 2 
I = I + l 
or TC 1 

2 TERM = PROC»A 
SUM = 0.0 

DO 3 M=1,N 

TERM = TEPM*(K-M+1 )*C/! ! J-M+1 )«B) 

3 SUM = SUM+TERM 
SIGF = SUM 
RFTUPN 

END 



FUNCTION FI !A,e,L,M) 

IMPLICIT PEAL^8!A-H»0-Z) 

NC'T=6 

C1=DGAMMA!0.5D 00) 

C6=CFL0AT! L+M+1 ) 

C2=DGAMMA!C6) 

C3=DGAMMA!L+1.5C 00) 

CA=2.0**! L+1 ) 

E=0.5*!L+M+1) 

C5=B**L/! A**2+B**2)**E 
F=0. 5*! L+l-M) 

G=0. 5«! 2*L+3) 

Z=B**2/ ( A**2+B**2 ) 

SUM=1.0 
TEPM=1.0 
OC 1 1=1,70 

TEPM=TERM*!E+I-1 ) *!F+I-1 )«Z/! !G+I-1 )♦! ) 
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IF(DABS(TEPM/SUM).LT.1.0D-09) GO TO 3 

1 SUM=SUM+TFPM 

2 WPITF(N0T,102) 

102 FCPMAT(T26t ' 70 TERMS') 

3 F1=C1^C2*C5/(C3*C4)*SUM 
RETURN 

END 



SUBROUTINE RI{JtM,K»KAPPA,RIl,RI2»P) 

IMPLICIT RFAL*8( A-H,0-Z) 

REAL*8 K, KAPPA 
DIMENSION F2I25) ,F3(25) 

A1=FACTPL( J+2) 

A2=FACTRL(^'+2) 

C1=KAPPA+1 

C2=2*KAPPA+1 

Ml=M+3 

M2=M+2 

Jl=J+3 

J2=J+2 

IF( J-M) 11 ,11 t9 
9 DO 10 N=1.J1 

F2(N)=F1 (C2 ,K,0, J+M+A-N) 

10 F3(N)=Fl(C2fK,l,J+M+4-N) 

GO TO 13 

11 DO 12 N=l,M.l 
F2(N)=F1(C2,K,0, J+M+4-N) 

12 F3(N)=Fl(C2fK,l, J+M+4-K) 

13 SUM=0.0 
TERM=1.0/M1 
DO 1 N=1,M1 
TERM=TERM*(M1-N+1 )/Cl 

1 SUM=SUM+TERM*F2( N) 

TERM=1.0/J1 

DC 2 M=1,J1 
TEPM=TERM*( Jl-N+1 )/Cl 

2 SUM=SUM+P*TERM*F2 { N) 

TEPM=1.0/M2 

DO 3 N=l, M2 
TEPM=TFRM*(M2-N-»-1 )/Cl 

3 SUM=SUM-TEPM*F2(N) 

TERM=1.0/J2 

DC 4 N=1,J2 
TERM=TERM*( J2-N+1 ) /Cl 

4 SUM=SUM-P*TERM*F2 (N) 

PI l=8*SUM/{ 2. 0/KAPPA )*4{J+M) 

SUM=C.0 
TEPM=1.C/M1 
DO 5 N=1,M1 
T£RM=TERM*( Ml-N+1 )/Cl 

5 SUM=SUM+TEPM#F3{N) 

TEPM=1.0/J1 

DO 6 N=1,J1 
TEPM=TERM=M Jl-N+1 )/Cl 

6 SUM=SUM+P*TERM*F3 (N) 

TERM=1.0/M2 

DC 7 N=1,M2 
TERM=TERM*(H2-N+1 )/Cl 

7 SUM=SUM-TERM4F3(N) 

TFPM=1.0/J2 

DO 8 N=1,J2 
TEP.M=TERM*( J2-N+D/C1 

8 SUM=SUM-P*TERM*F3(N) 

RI2=-8*(SUM+A2*F1 ( KAPPA, K , 1 , J ) /C1**M 1 +P*A 1*F 1 ( KAP PA,K , 
11 ,M) /C14+J1 ) /(2.0/KAPPA)**( J+M) 

RETURN 

END 



SUBROUTINE MIJ ( J , H , F.C.K , KAPPA t MTX IJ , S ) 
IMPLICIT REAL>«‘8(A-H,0-i) 
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INTEGER PfO 

PEAL^e K, KAPPA, MTXIJ 

A1=FACTRL( J+P+2) 

A2*FACTRL(M+Q+2) 

A3=FACTRL( J+0+2) 

AA=FACTRL(N+P+2) 

Cl= 2*KAPFA 
C2= 4*KAPFA 
Nl= N+O+J+F+5 
N2= N+O+J+P+4 
N3= M+O+J+P+6 
Fl= C1«^N1 
F2= C2**N1 
F3= C1**N2 
F4= C1**N3 

Sl= SIGF(P+Q+3,A2,N2,P+Q+3,C1,C2)/F2 
S2= SIGFU + P + 3 ,A1 ,N2, J+P+3,C1,C2 I/F2 
S3= SIGF(M+P+3 ,A4,N2,M+P + 3,C1,C2I/F2 
S4= SIGF( J + Q + 3 ,A3,N2 , J+C + 3 ,C1 ,C2 ) /F2 
Tl= (2*(K*>t'2-l )+A*KAPPA*<‘2)*(Al*A2+S*A3*A4)/F4 
T2=( -4*KAPFA*( P+1 ) )♦( A1*A2/ (J+P +2 ) +S* A3* A^/ ( M+P+2 ) ) / F 1 
T3=(-4*KAPPA*( 0+1 ) )*(Al*A2/(M+Q+2)+S*A3#A^/( J+0+2) )/Fl 
T4= 2*P*( P+l )*(Al*A2/( ( J+P+2 )*( J+P+1) )+S*A3*A4/( (M+p+2 
1)*(M+P+1) ) )/F3 

T5= 2*0*(C+1 )*( Al*A2 /( (M+q+2)*( P+0+1 ) ) +S*A3*A4/( (J+0+2 
1 »*( J+0+1) ) ) /F3 

MTXIJ=2*( T1+T2+T3+T4+T5+4*S1+4*S2+4*S*S3+4*S*S4)/ (2.0/ 
1KAPPA)**( J+M+P+Q) 

RFTURN 

END 



SUPRCUTINE 3I(K,B1,F) 

INPLICIT REAL*et A-H,0-Z) 

REALMS K 

NCT = e 

A1=K**2 

A2=1+A1 

A3=1-A1 

A4=A2**2 

A5=A2**3 

Tl=4/A‘5 

T2=4*A3/(A1*A4) 

T3=6*A3/(A1*A5) 

T4=2/(A1*A2) 

TERM=1.0 

SUM=TERM 

DC 1 N=2,7C 

Fl=( N*#2-1.0» /N**2 

TERN=-1*TEPM*F1*A1 

IF(DABS(TEPM/SUM),LT.1.0D-09) GC TO 3 

1 SUN=SUM+TERM 

2 WRITE(N0T,102) 

3 P1=SUN-P*T1+P*T2-P*T3+P*T4 

102 FCRNAT(T26, 'B1 DOES NOT CONVERGE IN 70 TERM'S') 
RETURN 
END 



SUBROUTINE BII(K,B2,P) 

IMPLICIT REAL*8(A-H,0-Z) 

REAL#8 K 

N0T=6 

A1=K**2 

A2=1+A1 

A3=1-A1 

A4=1+3*A1 

A5=K**3 

A6=A2**2 

A7=A2**3 

Tl=l/K 

T2=8*A3/(K*A7) 
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T3=e*K/A7 

Ta=12*A1/(A5»A6J 

T5=4*A4/(A5*A7) 

T6=4/(A5*A2) 

T7=4/( K+A6) 

TPRMl=5*K/6 

TERM2=^*A3/(K*A6» 

TFRM3=16*K/(3*A6) 

SHM1=TEPM1 
nr 1 N=3,70 

F2=( N-l,0)**2»(2.0*N+1.0)*(2.0*N-3.0)/( ( 2.0*N-1.0 )**2* 
1( N-2.0)*N) 

TERM1=-1*TERM1*F2*A1 

IF(nABS(TERMl/SUMl).LT.l,OD-09)GO TO 3 

1 SUM1=SUM1+TERM1 

2 WRITE(NOT,102» 

3 SUM=TERM2+TERM3 
DO 4 N=2t7C 
F3=(N-1.0)/N 

F4=(2.C*N-1,0) /(2.0*N+1.0) 

TERM2=-1*TERM2*F3=»A1 
TERM3=-1*TERM3*F4*A1 
TFRM=TERM2 + TER>'3 

IF(DABS(TEPM/SUM».LT.1.0D-09) GO TO 6 

4 SUM=SUM+TEPM 

5 WRITE(NOT» 102» 

6 B2=-l*( SUM1+P*SUM+T1-P*T2-P*T3+P*T4-P*T5 + P<‘T6-P*T7) 

102 F0RMAT(T26, *82 DOES NOT CONVERGE IN 70 TERMS') 

RETURN 

END 



SUBROUTINE BIII(KtB3,P) 

IMPLICIT REAL*8( A-H,0-Z) 

REAL*P K 

N0T=6 

A1=K»*2 

A2=1+A1 

A3=1+3*A1 

A4=2*A1 

A5=2+14*A1 

A6=A2**2 

A7=A2**3 

T1=3.141593/(3*K) 

T2=4/A7 

T3=6/(A1*A2) 

T4=A5/( A1*A7) 

T5=B/(A1*A6) 

T6=8/A6 
TERMl=Al/4 
SUM1=TERM1 
DO 1 N=3,70 

F5=(N-1.0)**2/( (N-2.0)*N) 

TEPM1=-1#TERM1*F5*A1 

IF(DABS(TEPMl/SUMl).LT.1.0D-09) GO TO 3 

1 SUM1=SUM1+TERM1 

2 WRITE(NOT,102) 

3 TEPM3=4*A3/( A1*A6) 

TERM4=8*A4/(3*A6) 

SUM=TERM3+TERM4 

DC 7 N=2.70 
F3=(N-1.0) /N 

F4=( 2.0*N-1.0) /(2.0+N+1.0) 

TERM3=-1*TERM3*F3*A1 

TEPM4=-1*TERM4*F4*A1 

TERM=TERM3+TERM4 

IF(0ABS(TEPM/SUM).LT.1.0D-09)G0 TO 9 
7 SUM=SUM+TERM 
a WRITE(N0T,102) 

9 B3=SUM1+P*SUM+T1+P*T2+P*T3-P*T4-P*T5-P*T6 
102 F0PMAT(T26, »B3 DOES NOT CONVERGE IN 70 TERMS') 
RETURN 
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END 



SUBROUTINE GAUSS 1 ( N , E PS , A , X , KER , K ) 

REAL*8 A,X,Y,D 

DIMENSION A( I ) ,X( 1 ) ,L( 100 ) ,M(100 ) ,Y ( 100, 100) 
DO 1 1=1, N 
DO 1 J=1,N 
IND=( I-l )*K+J 

1 Y(I,J)=A{IND) 

KEP = 1 
N2=2*N 

CALL AP.RAY1(2,N,N,100,100,Y,Y) 

CALL MINV1{Y,N,D,L,M) 

CALL ARRAY1(1,N,N,100,100,Y,Y) 

IFID.EO.O.) KER=2 
DO 2 1=1, N 
DO 2 J=1,N 
IND=( I-l KK + J 

2 X( IND)=Y(I , J) 

RETURN 

END 



SUBROUTINE ARR AYl ( MODE , I , J,N,M,S,D) 
DIMENSION S(l) ,D(1) 

REAL*8 S,D 
NI=N-I 

IF(MCDE-l) 100,100,120 
100 IJ=I*J+1 
NM=N*J+1 
DC 110 K=1,J 
NM=NM-MI 
DO 110 L=1 ,I 
I J=I J-1 
NM=NM-1 

no 0(NM)=S(IJ) 

GO TO lAO 
12C IJ=0 
NM=0 

DC 130 K = 1 , J 
DO 125 L=1,I 
U = I J+1 

NM=NM+1 

125 S(IJ)=D(NM) 

130 NM=NM+NI 
lAO RETURN 
END 



SUBROUTINE MINVl ( A,N,C,L,M) 
DIMENSION A(1),L(1),M(1) 

DOUBLE PRECISION A , D , B IGA ,HOLD 

D=1.0 

NK=-N 

DC 80 K=1,N 
mk=nk+n 

L(K)=K 

M(K)=K 

KK=NK+K 

BIGA=A(KK) 

DO 20 J=K,N 
I Z=N*( J-1 ) 

DO 20 I=K,N 
IJ=IZ+I 

10 IF(DABS(BIGA)-DABS( A( IJ) ) ) 15,20,20 
15 BIGA=A(IJ) 

L(K)=I 
M(K)=J 
20 CONTINUE 
J=L( K) 

IF(J-K) 35,35,25 
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25 KI=K-N 

DO 30 I=1,N 
KI=KI+N 
HOLD=-A( KI ) 

JI=KI-K+J 
A(KI ) = A( JI ) 

3C A( JI »=HOLD 
35 I=M(K) 

IF(I-K) 45,45»38 
38 JF=N*(I-1) 

DC AO J=1»N 
JK=NK+J 
JI=JP+J 
HCLD=-A( JK) 

A( JK)=A( JI ) 

40 A( JI »=HOLD 

45 IF(BIGA) 48,46,48 

46 D=0.0 
RETURN 

48 DC 55 1=1, N 

IF(I-K) 50,55,50 
50 IK=NK+I 

A( IK)=A( IK) /(-BIGA) 

55 CONTINUE 
DO 65 1=1, N 
IK=NK+I 
I J=I-N 
DO 65 J=1,N 
I J=I J+N 

ini-K) 60,65,60 
60 IF(J-K) 62,65,62 
62 KJ=IJ-I+K 

A(IJ)=A( IK )*A(KJ)+A(I J) 
65 CONTINUE 
KJ=K-N 
DC 75 J=1,N 
KJ=K J+N 

IF(J-K) 70,75,70 
70 A(KJ)=A(KJ)/8IGA 
75 CONTINUE 

A( KK ) =1.0/BIGA 
80 CONTINUE 
K=N 

100 K=K-1 

IF(K) 150,150,105 
105 I=L(K) 

IF(I-K) 120,120,108 
108 J0=N*(K-1) 

JR=N*( T-1) 

DC 110 J=1,N 
JK=JO+J 
HCLD=A( JK) 

JI=JR+J 
A( JK)=-A( JI ) 

110 A( JI )=HOLD 
120 J=M(K) 

IF(J-K) 100,100,125 
125 KI=K-N 

DO 130 1=1, N 
KI=KI+N 
HCLD=A( KI ) 

JI=KI-K+J 
A( KI ) =-A( JI ) 

130 A( JI )=HOLD 
GC TO 100 
150 RETURN 
END 
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